#usage: python qe_exp_time_cal.py op tabdir 400 50000 60 outdir
import numpy as np
from astropy.table import Table
import sys

band = sys.argv[1] #uv op or ir
tabdir = sys.argv[2] #theoretical QE table directory
wl = float(sys.argv[3]) # measured wavelength
dn = float(sys.argv[4]) # measured DN
bias = float(sys.argv[5]) # measured bias
exptime = float(sys.argv[6]) # exposure time for measured wavelength
outdir = sys.argv[7] # output directory

wave_ratio = np.array([280,320,360,400,440,480,520,560,600,640,680,720,760,800,826.87,862.21,887.30,902.31,908.02,914.06,920.98,935.69,946.02,965.72,975.16,991.08,\
997.81,1001.80,1040.12,1080,1081.39])
ratio = np.array([23.41,25.07,25.46,25.76,25.47,25.99,24.77,24.56,24.23,23.71,23.32,22.41,22.45,21.48,29.98,21.95,31.52,27.10,29.80,28.01,30.28,24.28,24.64,22.96,\
43.08,26.38,28.96,27.07,17.26,15.88,17.15])
dn -= bias
DNmax = 50000/1.3-bias  #max 50000e-
DNmin = 5000-bias  #min 5000 DN

if band == 'op':
    tab = Table.read(tabdir+"/qe_op.tab")
elif band == 'uv':
    tab = Table.read(tabdir+"/qe_uv.tab")
elif band == 'ir':
    tab = Table.read(tabdir+"/qe_ir.tab")
else:
    print("invalid waveband!")

wave_out = np.arange((230,1030,10))
wave_qe = tab['wave']
qe = tab['qe_corr']

qe_theory = np.interp(wave_out,wave_qe,qe)
qe_theory /= np.interp(wl,wave_qe,qe)
ratio_theory = np.interp(wave_out,wave_qe,qe)
ratio_theory /= np.interp(wl,wave_ratio,ratio)
exp_times_max =  (DNmax/dn)/(qe_theory*ratio_theory)*exptime
exp_times_min =  (DNmin/dn)/(qe_theory*ratio_theory)*exptime
result_tab = Table()
result_tab['wave'] = wave_out
result_tab['wave'].unit = 'nm'
result_tab['exp_time_min'] = exp_times_max
result_tab['exp_time_min'].unit = 's'
result_tab['exp_time_max'] = exp_times_min
result_tab['exp_time_max'].unit = 's'
result_tab.writeto(outdir+"/exp_times.tab",format='ipac',overwrite=True)
